# Script for processing West Coast Trajectories (2/2019)
# Takes raw tracks, breaks at ports and stores port information

#### This code can be crudely parallelized by running each year as a separate process (e.g. different command window)
# C:\Python27\ArcGISx6410.5\python.exe path_of_file\process_linesM_fuel_07032019.py

# for use with anaconda
try:
    import archook #The module which locates arcgis
    archook.get_arcpy()
    import arcpy
except ImportError:
    error
    # do whatever you do if arcpy isnt there.


#import arcpy
#import numpy
import datetime
import time
import os

#import fuel_calculator 

# run using: C:\Python27\ArcGISx6410.5\python.exe C:\Users\rklotz\OneDrive\Work\Ports_ECAs\arcpy_tools\process_linesM_fuel_07032019.py

arcpy.env.overwriteOutput=True
#numpy.set_printoptions(threshold=numpy.nan)

### California ECA shape file
CAeca2009_file = r"H:\My Drive\Boats\ReplicationCode\data\boundaries\CAeca_200907.shp"
CAeca2011_file = r"H:\My Drive\Boats\ReplicationCode\data\boundaries\CAeca_201112.shp"

portbound_file = r"H:\My Drive\Boats\ReplicationCode\data\port_buffers.gdb\west_coast_ports_eca100_AK_HI_v3" 
dis_field = "OBJECTID"
boundary_cond = 'PORT_ID>100 and PORT_ID<120' # port buffers that define boundary

## input and output options #############################
lines_file_tag = r"H:\My Drive\Boats\ReplicationCode\data\AIS\lines_v2\merged_lines.gdb\cleaned_v2_%s" # these include interpolated separately
out_gdb = r"H:\My Drive\Boats\ReplicationCode\data\AIS\voyages.gdb"

# fields for 2009-2014
flds_to_copy_0914 = ['MMSI', 'TrackStartTime', 'TrackEndTime', 'IMO', 'Name', 'VesselType', 'Length', 'Width', 'DimensionComponents','OUT_CAECA','OUT_NAECA','INTERP_FLAG','TRACK_ID','INTERP_TRACKS']

# use this one for 2015-16 because available data fields are different
flds_to_copy_1516 = ['MMSI', 'TrackStartTime', 'TrackEndTime','IMO' , 'VesselName', 'VesselType', 'Length', 'Width', 'Draft','OUT_CAECA','OUT_NAECA','INTERP_FLAG','TRACK_ID','INTERP_TRACKS']

#line_selector = 'not VesselType<60 and not VesselType=99'
# selecting any that are classified as cargo or tanker
# then getting large passenger ships
# then getting large other ships (either no vessel type or "other type")
line_selector_in = '( (VesselType Between 70 and 89 or VesselType In (1003, 1004, 1016, 1017, 1024))  or (( VesselType Between 60 and 69 or VesselType Between 1012 and 1015)  and Length>=60) or ((VesselType Between 90 and 999 or VesselType=0) and Length>=60 ) )'


# using GCS_WGS_1984
spatial_reference = arcpy.SpatialReference(4326) # GCS_WGS_1984

years = range(2009,2017)

# testing
years = [2009]

def main() :

    # Create output gdb if needed
    if not arcpy.Exists(out_gdb) == True :
        print "Create Output GDB"
        #print os.path.dirname(out_gdb) , os.path.basename(out_gdb)
        arcpy.CreateFileGDB_management( os.path.dirname(out_gdb) , os.path.basename(out_gdb) )
    
    # run code twice, once splitting at study area boundary and once not
    for split_at_boundary in [1,0] :
        if split_at_boundary==0:
            out_fc_tag = os.path.join(out_gdb,"wc_no_bound_%s") 
        else:
            out_fc_tag = os.path.join(out_gdb,"wc_%s")     
        print out_fc_tag
        
        for year in years :
    
            lines_file = lines_file_tag % (year) 
            out_fc = out_fc_tag % (year) 
    
            if year in [2009,2010,2011,2012,2013,2014]:
                flds_to_copy = flds_to_copy_0914
            else:
                flds_to_copy = flds_to_copy_1516
                
            #### getting lines feature class (selecting tracks)
    
            # dropping specific vessels
            if year>2009 and year<2015 :
                line_selector = line_selector_in + "and not MMSI In (366240119,366907020,366194506,368080090,366758000,636461104,366208968,312880300,367163760,366978210,367119810,367151450)"
                # 366907020 cargo ship travels circles really close to port boundary
                # 369080090 tanker with complicate geometry
                # 366972000 seems like a dredge
                # 366194506 is a cargo vessel that operates just off LA
                # 636461104 is cargo vessel with complicated geometry
                # 312880300 entrance with complicated geometry (fishing)
                # 367163760,366978210,367119810,367151450 ferries in gulf 
    
                # east coast line in august 2012 that stalls on intersect
                line_selector = line_selector + "and not MMSI in (366914210,369005505,366108927,367407028,366904983,367504802,367730661,366108927,367911018,367045511)"
                # 367407028 appears to be a ferry from Conn to LI
                # 366904983 ferry from Conn to LI
                # 367730661,366108927,367911018,367045511 ferries in gulf (galveston)
            
                # 366914210	366240119	M/V Deleware
    
    
            
            else: # 2009, 2015, 2016
                line_selector = line_selector_in + "and not MMSI In (316146000,338570000,366972000)"
                # 366972000 is a dredge
                
    
    ##        if year==2009 :
    ##            # canadian military tanker
    ##            line_selector = line_selector "and not MMSI In (316146000)"
    ##
    ##        if year==2010 :
    ##            line_selector = line_selector + "and not MMSI In (366907020, 366194506)" # 366907020 really complicate track (maybe research vessel)
    ##                                                                                     # 366194506 is a cargo vessel that operates just off LA
    ##        if year==2011 :
    ##            # tanker with very complicated entrances/exits
    ##            line_selector = line_selector "and not MMSI In (366907020)"
    ##
    ##        if year==2013 :
    ##            # cargo ship travels circles really close to port boundary
    ##            # tanker with complicate geometry
    ##            # seems like a dredge
    ##            line_selector = line_selector + "and not MMSI in (366907020,368080090, 366758000)"
    ##            
    ##        if year==2016 :
    ##            # appears to be dredging south of seattle 
    ##            line_selector = line_selector + "and not MMSI=338570000 and not MMSI=366972000"
    
            # do splitting steps
            processLines(lines_file,portbound_file,out_fc,line_selector,flds_to_copy)

    

def processLines(lines_file,portbound_file,out_fc,line_selector,flds_to_copy) :

    start = time.time()
    print "make feature layer"
    arcpy.management.MakeFeatureLayer(lines_file,"lines_fc",line_selector) 
    
    # make port feature layer
    # this is used to determine where endpoints of tracks are
    # need to remove boundary if not breaking there
    port_selector = 'not ( PORT_ID>=50 and PORT_ID<60 )' # ignoring any that are waypoints (in 50s)
    if split_at_boundary==0 :
        port_selector = port_selector + "and not (" + boundary_cond + ")" # removing boundaries
              
    arcpy.management.MakeFeatureLayer(portbound_file,"portbound_int_fc",port_selector) 

    split_dir = "in_memory" #r"E:/data/temp.gdb"
    split_name = r"splittemp"
    split_fc = os.path.join(split_dir,split_name)

    # create feature class for broken lines
    #arcpy.CreateFeatureclass_management(split_dir, split_name,geometry_type="POLYLINE",spatial_reference=spatial_reference,has_m = "ENABLED")
    arcpy.CreateFeatureclass_management(split_dir, split_name,geometry_type="POLYLINE",spatial_reference=spatial_reference,has_m = "ENABLED",template="lines_fc")

    arcpy.AddField_management(split_fc,"PORT1","SHORT")
    arcpy.AddField_management(split_fc,"PORT2","SHORT")
    arcpy.AddField_management(split_fc,"TIME1","DATE")
    arcpy.AddField_management(split_fc,"TIME2","DATE")
    arcpy.AddField_management(split_fc,"RIGHT1","SHORT")
    arcpy.AddField_management(split_fc,"RIGHT2","SHORT")

    arcpy.AddField_management(split_fc,"DIST","FLOAT")
    arcpy.AddField_management(split_fc,"DIST_ECA2009","FLOAT")
    arcpy.AddField_management(split_fc,"DIST_ECA2011","FLOAT")
    arcpy.AddField_management(split_fc,"TIME_ECA2009","FLOAT")
    arcpy.AddField_management(split_fc,"TIME_ECA2011","FLOAT")

    
    #flds = arcpy.ListFields(split_fc)
    #field_names = [f.name for f in flds]

    #flds_to_copy = ['MMSI', 'TrackStartTime', 'TrackEndTime', 'IMO', 'CallSign', 'Name', 'VesselType', 'Length', 'Width', 'DimensionComponents']

    # fields to add after splitting by port and ECA boundaries
    flds_port = ["PORT1","TIME1","RIGHT1","PORT2","TIME2","RIGHT2","DIST_ECA2009","TIME_ECA2009","DIST_ECA2011","TIME_ECA2011"]

    # converts feature to geometry
    #port_geom = arcpy.CopyFeatures_management("portbound_int_fc",arcpy.Geometry()) #[0]

    # get CA ECA geometries
    CAeca2009_geom = arcpy.CopyFeatures_management(CAeca2009_file,arcpy.Geometry())[0]
    CAeca2011_geom = arcpy.CopyFeatures_management(CAeca2011_file,arcpy.Geometry())[0]

    # create single study area boundary
    us_bound_fc = "in_memory/usbound"
    arcpy.management.MakeFeatureLayer(portbound_file,"USbound_fc",boundary_cond)
    arcpy.Dissolve_management("USbound_fc",us_bound_fc,dissolve_field=dis_field) 
    USbound_geom = arcpy.CopyFeatures_management(us_bound_fc,arcpy.Geometry())
    arcpy.DeleteFeatures_management(us_bound_fc)
    
    # check to make sure USbound is a single line
    if len(USbound_geom)==1 :
        USbound_geom = USbound_geom[0]
    else :
        problem_with_USbound # cheap exception
            
    # getting port geometry and port ids
    port_geom = []
    port_ids = []
    port_geom_port = []
    port_ids_port = []
    port_geom_US = []
    port_ids_US = [] 
    for row in arcpy.da.SearchCursor("portbound_int_fc", ["SHAPE@","PORT_ID"],spatial_reference=spatial_reference) :
        port_geom.append(row[0])
        port_ids.append(row[1])
        if row[1]<100 :
            port_geom_port.append(row[0])
            port_ids_port.append(row[1])
        else :
            port_geom_US.append(row[0])
            port_ids_US.append(row[1])
            
    del row

    # get indices of port ids<100 and intersecting points
    ##                intind_port = [i for i in intind if port_ids[i]<100]
    ##                intpoints_port = [intpoints[i] for i in range(0,len(intind)) if port_ids[intind[i]]<100]

    #PORTS = len(port_ids)

    cur = None
    try:
        # Open an insert cursor for the new feature class
        cur = arcpy.da.InsertCursor(split_fc, ['SHAPE@',"DIST"] + flds_to_copy + flds_port )

        #### looping through lines and breaking
        print "looping through lines and breaking"
        splitlines = []
        skipped_complicated = 0 
        for row in arcpy.da.SearchCursor("lines_fc", ['SHAPE@'] + flds_to_copy ,spatial_reference=spatial_reference ) :
            line = row[0]

            #print("Feature {}:".format(row[2]))
            print(u'{0}, {1}'.format(row[1],row[2]))
            #print line.pointCount

            # only doing lines
            if line is None :
                print "skipped"
            else:
            #if line is not None :

                # drop condition based on geometry
                # dropping if length of segment is 0 or if geometry is very dense (10 points per km of length)
                #num_points = line.pointCount
                length = line.getLength('GEODESIC','KILOMETERS')

                skip_geom = 0 
                if length==0 :
                    skip_geom = 1
                else :
                    pt_density = line.pointCount/line.getLength('GEODESIC','KILOMETERS')
                    if pt_density>10 :
                        skip_geom = 1
                    print pt_density

                if skip_geom==1 :
                    print "skipped due to complicated geometry"
                    skipped_complicated+=1
                else:
                
                    # finding points of intersection with ports   
                    #print "intersect with ports"
                    intpoints_port, intind_port = get_int_points(line,port_geom_port)

                    #print intpoints_port

                    if len(intpoints_port)>0 :
                        # lines that intersect ports
                        
                        # split at ports
                        #print "split at ports"
                        splitline=arcpy.SplitLineAtPoint_management(line, intpoints_port, arcpy.Geometry(),"1 meters")

                        # loop through split lines
                        for ln in splitline:
                            # get port ids and right/left for each lines endpoints
                            port_list_ln = get_start_end(ln,port_geom_port,port_ids_port)

                            if split_at_boundary==1 :

                                # test if endpoints left of boundary
                                ptfirst = arcpy.PointGeometry(ln.firstPoint,spatial_reference)
                                ptlast = arcpy.PointGeometry(ln.lastPoint,spatial_reference)
    
                                qpd_first = USbound_geom.queryPointAndDistance(ptfirst)
                                qpd_last = USbound_geom.queryPointAndDistance(ptlast)
    
                                # determine if line terminates outside US boundary
                                # if either endpoint is not classified AND falls to left of US bound (False)
                                if ( port_list_ln[0] is None and qpd_first[3] is False ) or ( port_list_ln[3] is None and qpd_last[3] is False ):
                                    # if so, try to split at US boundary
    
                                    # finding points of intersection with US boundary   
                                    intpoints_US, intind_US = get_int_points(ln,port_geom_US)
                        
                                    if len(intpoints_US)>0 :
                                        #print "split at boundaries
                                        splitline2=arcpy.SplitLineAtPoint_management(ln, intpoints_US, arcpy.Geometry(),"1 meters")
    
                                        # add lines split at US boundary
                                        for ln2 in splitline2 :
                                            port_list_ln2 = get_start_end(ln2,port_geom,port_ids)
                                            eca_list = get_eca_data(ln2,CAeca2009_geom,CAeca2011_geom)
                                            ##### ADDING
                                            cur.insertRow( [ln2,ln2.getLength('GEODESIC','KILOMETERS')]+list(row[1:]) + port_list_ln2 + eca_list )
                                    else :
                                        # if no intersection, add unsplit line
                                        ##### ADDING
                                        eca_list = get_eca_data(ln,CAeca2009_geom,CAeca2011_geom)
                                        cur.insertRow( [ln,ln.getLength('GEODESIC','KILOMETERS')]+list(row[1:]) + port_list_ln + eca_list )
                                    
                                else :
                                    # if line doesn't terminate outside, add segement only split by port to table
                                    ##### ADDING
                                    eca_list = get_eca_data(ln,CAeca2009_geom,CAeca2011_geom)
                                    cur.insertRow( [ln,ln.getLength('GEODESIC','KILOMETERS')]+list(row[1:]) + port_list_ln + eca_list )
                
                            else :
                                # not splitting at boundary 
                                eca_list = get_eca_data(ln,CAeca2009_geom,CAeca2011_geom)
                                cur.insertRow( [ln,ln.getLength('GEODESIC','KILOMETERS')]+list(row[1:]) + port_list_ln + eca_list )
                    else :
                        # if lines don't intersect ports

                        if split_at_boundary==1 :
                            # finding points of intersection with US boundary   
                            intpoints_US, intind_US = get_int_points(line,port_geom_US)
                                    
                            # if intersection split at US boundary
                            if len(intpoints_US)>0 :
                                splitline2=arcpy.SplitLineAtPoint_management(line, intpoints_US, arcpy.Geometry(),"1 meters")
    
                                # add lines split at US boundary
                                for line2 in splitline2 :
                                    port_list_line2 = get_start_end(line2,port_geom,port_ids)
                                    ##### ADDING
                                    eca_list = get_eca_data(line2,CAeca2009_geom,CAeca2011_geom)
                                    cur.insertRow( [line2,line2.getLength('GEODESIC','KILOMETERS')]+list(row[1:]) + port_list_line2 + eca_list )
                        else :
                            # if not breaking at study area boundary 
                            eca_list = get_eca_data(line,CAeca2009_geom,CAeca2011_geom)
                            port_list_line = get_start_end(line,port_geom,port_ids)
                            cur.insertRow( [line,line.getLength('GEODESIC','KILOMETERS')]+list(row[1:]) + port_list_line + eca_list )
                            


    except Exception as e:
        print "failed in loop through lines"
        print(e)
    finally:
        # Cleanup the cursor if necessary
        if cur:
            del cur
            
    print("--- %s seconds ---" % (time.time() - start))
    print(u'{0}'.format(skipped_complicated))
    
    # add id for split tracks - calling it voyage
    arcpy.AddField_management(split_fc,"VOYAGE_ID","LONG")
    arcpy.CalculateField_management(split_fc, "VOYAGE_ID","!OBJECTID!","PYTHON")
    
    #### intersecting lines with port buffers       
    fms = arcpy.FieldMappings() # creating field mappings
    fms.addTable(split_fc)      # adding all fields from lines

    fm1 = arcpy.FieldMap()          # adding field map for getting port IDs
    fm1.addInputField(portbound_file,'PORT_ID')

    # creating outfield
    outfield = arcpy.Field()
    outfield.name = "PORT_ID"
    outfield.type = "String"
    outfield.length = 100
    outfield.scale = 0

    fm1.outputField = outfield
    fm1.mergeRule = "Join"
    fm1.joinDelimiter = ","

    fms.addFieldMap(fm1) # adding additional outfield

    #match_option = "INTERSECT"
    #search_radius = "2 meters"

    match_option = "WITHIN_A_DISTANCE"
    search_radius = "5 kilometers"

    print "joining port buffers"
    # Replace a layer/table view name with a path to a dataset (which can be a layer file) or create the layer/table view within the script
    # The following inputs are layers or table views: "westcoast_lines_2009_1", "west_coast_ports_v6"
    arcpy.SpatialJoin_analysis(target_features=split_fc, join_features=portbound_file, out_feature_class=out_fc
                , join_operation="JOIN_ONE_TO_ONE", join_type="KEEP_ALL", field_mapping=fms
                , match_option=match_option, search_radius=search_radius, distance_field_name="")

    print("--- %s seconds ---" % (time.time() - start))

    # cleaning up
    #arcpy.DeleteFeatures_management(intpoints_fc)
    arcpy.DeleteFeatures_management(split_fc)
    

def eca_intersect(ln,eca) :
    dist = 0 
    time = 0
    
    # check line to see if not disjoint with ECA boundary
    if not ln.disjoint(eca) :
        #ln_eca = ln.intersect(eca, 2) # doesn't work

        # clipping (get array of a single multipart line)
        ln_eca_array = arcpy.Clip_analysis(ln,eca,arcpy.Geometry())

        #print ln_eca_array
        #asdasd

        # if array is not empty
        if len(ln_eca_array)>0 :
            ln_eca = ln_eca_array[0]
               
            # calculate distance in ECA
            dist = ln_eca.getLength('GEODESIC','KILOMETERS')
            #print dist

            # calculate time in ECA
            # loop through all parts of the multipart line
            # calculate time between first and last point
            # sum

            time = 0 
            for i in range(ln_eca.partCount) :
                prt = ln_eca.getPart(i)

                tfirst = datetime.datetime.fromtimestamp(prt[0].M) 
                tlast = datetime.datetime.fromtimestamp(prt[-1].M)

        ##            print i 
        ##            print prt[0]
        ##            print prt[-1]
        ##
        ##            print tfirst
        ##            print tlast
                
                
                dt = tlast - tfirst

                # summing time across all intervals
                # (converting seconds to hours)   
                time += abs( dt.total_seconds() / 3600 )

                #print time
    
    return dist, time


def get_eca_data(ln,CAeca2009_geom,CAeca2011_geom) :

    DIST_ECA2009,TIME_ECA2009 = eca_intersect(ln,CAeca2009_geom)
    DIST_ECA2011,TIME_ECA2011 = eca_intersect(ln,CAeca2011_geom)

    #create list    
    eca_int_list = [DIST_ECA2009,TIME_ECA2009,DIST_ECA2011,TIME_ECA2011]
   
    return eca_int_list

def get_int_points(line,port_geom): 
    # finding points of intersection (looping through all points)
    
    PORTS = len(port_geom)
    
    intpoints = []
    intind = [] # index of ports that intersected
    for port_i in range(0,PORTS) :
        port_line = port_geom[port_i]
        #print port_i
        if not line.disjoint(port_line) :
            intersect = line.intersect(port_line, 1)
            intpoints.append(intersect)
            intind.append(port_i)

    return intpoints, intind

def get_start_end(ln,port_geom,port_ids) :
    PORTS = len(port_ids)

    # get date/time for first and last points
    tfirst = datetime.datetime.fromtimestamp(ln.firstPoint.M) 
    tlast = datetime.datetime.fromtimestamp(ln.lastPoint.M)

    ptfirst = arcpy.PointGeometry(ln.firstPoint,spatial_reference)
    ptlast = arcpy.PointGeometry(ln.lastPoint,spatial_reference)

    # get point 
    ptfirst2 = arcpy.PointGeometry(ln.getPart(0)[1],spatial_reference)
    ptlast2 = arcpy.PointGeometry(ln.getPart(0)[-2],spatial_reference)

    # looping through all ports
    # finding port that is closest to first and last point
    dist_first = 1e9
    dist_last = 1e9
    id_first = None
    id_last = None
    right_first = None
    right_last = None 
    for port_i in range(0,PORTS) : # portint_ind :
        #print port_ids[port_i]
        
        # getting distance between first and last
        dist_first_tmp = ptfirst.distanceTo(port_geom[port_i])
        dist_last_tmp = ptlast.distanceTo(port_geom[port_i])
        
        # Using this to get right vs left (second, and second last points)
        qpd_first2 = port_geom[port_i].queryPointAndDistance(ptfirst2)
        qpd_last2 = port_geom[port_i].queryPointAndDistance(ptlast2)

        #dist_first_tmp = qpd_first[2]
        #dist_last_tmp = qpd_last[2]
        
        # update
        if dist_first_tmp < dist_first :
            id_first = port_ids[port_i]
            dist_first = dist_first_tmp
            right_first = int( qpd_first2[3] )

        # update
        if dist_last_tmp < dist_last :
            id_last = port_ids[port_i]
            dist_last = dist_last_tmp
            right_last = int( qpd_last2[3] )

    # 1e-3 should be within about 100 meters
    # 3e-2 should be within 3km
    dist_buff = 3e-2
    if dist_first>dist_buff :
        id_first = None
        right_first = None
        
    if dist_last>dist_buff :
        id_last = None
        right_last = None 
        

    if tlast<tfirst :
        # last point occurred first in terms of time
        port_list = [ id_last , tlast , right_last , id_first, tfirst, right_first ]
    else :
        # first point occurred first
        port_list = [ id_first, tfirst , right_first , id_last , tlast , right_last]

    return port_list

if __name__ == "__main__" :
    main() 



#### THIS WORKS FOR BREAKING LINES BUT IS RESOURCE INTENSIVE
###### break lines at port buffers
##print "get points of intersection"
###intpoints_fc = r"in_memory\intpoints"
##intpoints_fc = r"E:\data\temp.gdb\intpoints"
##arcpy.Intersect_analysis([ "portbound_int_fc" ,"lines_fc" ], intpoints_fc, output_type="POINT")
##print("--- %s seconds ---" % (time.time() - start))
##
###print "Integrating"
###arcpy.Integrate_management (["lines_fc",intpoints_fc], "0.1 feet")
##
##print "break lines at points"
###broke_lines_fc = r"in_memory\broke_lines"
##broke_lines_fc = r"E:\data\temp.gdb\broke_lines"
##arcpy.SplitLineAtPoint_management("lines_fc", intpoints_fc, broke_lines_fc,"1 meters")
##print("--- %s seconds ---" % (time.time() - start))
